Light-cone-like spreading of correlations in a quantum many-body system 
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How fast can correlations spread in a quantum many-body system? Based on the seminal work 
by Lieb and Robinson [T], it has recently been shown that several interacting many-body systems 
exhibit an effective light cone that bounds the propagation speed of correlations [2 - 5 . The existence 
of such a "speed of light" has profound implications for condensed matter physics and quantum 
information, but has never been observed experimentally. Here we report on the time-resolved 
detection of propagating correlations in an interacting quantum many-body system. By quenching 
a one-dimensional quantum gas in an optical lattice, we reveal how quasiparticle pairs transport 
correlations with a finite velocity across the system, resulting in an effective light cone for the 
quantum dynamics. Our results open important perspectives for understanding relaxation of closed 
quantum systems far from equilibrium [6] as well as for engineering efficient quantum channels 
necessary for fast quantum computations [7]. 



In contrast to relativistic quantum field theory, no 
"speed limit" exists in non-relativistic quantum mechan- 
ics, allowing in principle for the propagation of informa- 
tion over arbitrary distances in arbitrary short times [2]. 
However, one could naively expect that in real physi- 
cal systems short-range interactions allow information to 
propagate only with a finite velocity. The existence of a 
maximal velocity, also called Lieb-Robinson bound, has 
indeed been shown theoretically in some systems, e.g. in- 
teracting spins on a lattice [2H5], but to which extent this 
result can be generalised remains an open question [8HTT]. 
Lieb-Robinson bounds have already found a number of 
fundamental applications [12) H3]. For example, they al- 
low for a rigorous proof of a long-standing conjecture that 
linked the presence of a spectral gap in a lattice system 
to the exponential decay of correlations in the ground 
state [HJ [15]. They also provide fundamental scaling 
laws for the entanglement entropy, which is an indicator 
of the computational cost for simulating strongly inter- 
acting systems [T5j . 

In the context of quantum many-body systems, the 
existence of a Lieb-Robinson bound can be probed by 
recording the dynamics following a sudden parameter 
change (quench) in the Hamiltonian. In that case, a 
simple picture has been suggested: quantum-entangled 
quasiparticles emerge from the initially highly excited 
state and propagate ballistically [3j, carrying correla- 
tions across the system. Ultracold atomic gases offer 
an ideal testbed to explore such quantum dynamics due 
to their almost perfect decoupling from the environment 
and their fast tunability [17]. In addition, the recently 
demonstrated technique of single-site imaging in an op- 
tical lattice [T8l [19] offers the resolution and sensitivity 
necessary to reveal the dynamical evolution of a many- 
body system at the single-particle level. 



Our system consists of ultracold bosonic atoms in an 
optical lattice and is well described by the Bose-Hubbard 
model [201 |2T]. This model is parameterised by two en- 
ergy scales: the on-site interaction, U, and the tunnel 
coupling between adjacent sites, J. Driven by the com- 
petition of these two parameters, a quantum phase tran- 
sition between a superfluid and a Mott-insulating phase 
occurs in homogeneous systems with integer filling n. In 
the one-dimensional (Id) geometry considered here, the 
critical point of this transition is located at (U/J) c ~ 3.4 
[22] . We observed the time evolution of spatial corre- 
lations after a fast decrease of the effective interaction 
strength U/J from an initial value deep in the Mott- 
insulating regime, with filling n = 1, to a final value closer 
to the critical point (Fig. [IJi). After such a quench, the 
initial many-body state |^o) is highly excited and acts 
as a source of quasiparticles. In order to elucidate the 
nature and the dynamics of these quasiparticles, we have 
developed an analytical model in which the occupancy 
of each lattice site is restricted to n = 0, 1 or 2 (see 
Appendix). For large interaction strengths, the quasi- 
particles consist of either an excess particle (doublon) 
or a hole (holon) on top of the unity-filling background. 
Fermionizing these quasiparticles with a Jordan-Wigner 
transformation allows us to partially eliminate the non- 
physical states in which a lattice site would be occupied 
by two quasiparticles. To first order in J/U, we then 
find that the many-body state at time t after the quench 
reads: 

|¥(t)> ~ |* >+ iV8^J2 sin ( ka ^ 



*Mfc)+«h(-fc)]t/fil $ k hl k \* ) , (1) 



with ai a t the lattice period. Here d\ and hi are the 



2 



* A A A A A A A (r 



quench 




position 



\%l \#/ w w \m/ \m/ \9/ \m/ 



d = vt 



FIG. 1. Spreading of correlations in a quenched atomic 
Mott insulator, a, A Id ultracold gas of bosonic atoms 
(black balls) in an optical lattice is initially prepared deep 
in the Mott-insulating phase with unity filling. The lattice 
depth is then abruptly lowered, bringing the system out of 
equilibrium, b, Following the quench, entangled quasiparticle 
pairs emerge at all sites. Each of these pairs consists of a 
doublon (red ball) and a holon (blue ball) on top of the unity- 
filling background, which propagate ballistically in opposite 
directions. It follows that a correlation in the parity of the 
site occupancy builds up at time t between any pair of sites 
separated by a distance d = vt, where v is the relative velocity 
of the doublons and holons. 



creation operators for a doublon and a holon with mo- 
mentum fc, respectively, and k belongs to the first Bril- 
louin zone. Quasiparticles thus emerge at any site in the 
form of entangled pairs, consisting of a doublon and a 
holon with opposite momenta. Some of these pairs are 
bound on nearest-neighbour sites while the others form 
wave packets, due to their peaked momentum distribu- 
tion. The wave packets propagate in opposite directions 
with a relative group velocity v determined by the dis- 
persion relation €d(k) + 6h( — k) of doublons and holons 
(Fig. [TJ)). The propagation of quasiparticle pairs is re- 
flected in the two-point parity correlation functions [23J: 

Cd(t) = (Sj^Sj+dit)) - (8j(t))(8 j+d (t)) , (2) 

where j labels the lattice sites. The operator sj(t) = 
e w[nj(t)-n] measures the parity of the occupation number 
fij{t). It yields +1 in the absence of quasiparticles (odd 
occupancy) and -1 if a quasiparticle is present (even occu- 
pancy). Because the initial state is close to a Fock state 
with one atom per lattice site, we expect Cd{t = 0) ~ 0. 
After the quench, the propagation of quasiparticle pairs 
with the relative velocity v results in a positive correla- 
tion between any pair of sites separated by a distance 
d = vt. 

The experimental sequence started with the prepara- 
tion of a two-dimensional (2d) degenerate gas of 8T Rb 
confined in a single antinode of a vertical optical lattice 
fT9l [23] (z-axis, a\ at = 532 nm). The system was then 



divided into about 10 decoupled Id chains by adding a 
second optical lattice along the ?/-axis and by setting both 
lattice depths to 20.0(5) E r , where E r = (2nh) 2 / (8mtf at ) 
is the recoil energy of the lattice and m the atomic mass 
of 87 Rb. The effective interaction strength along the 
chains was tuned via a third optical lattice along the 
x-axis. The number of atoms per chain ranged between 
10 and 18, resulting in a lattice filling ft = 1 in the Mott- 
insulating domain. The inital state was prepared by adi- 
abatically increasing the x-lattice depth until the interac- 
tion strength reached a value of (U/J)o = 40(2). At this 
point, we measured the temperature to be T ~ 0.1 U/k B 
(k B is the Boltzmann constant) following the method de- 
scribed in Ref. [19]. We then brought the system out of 
equilibrium by lowering the lattice depth typically within 
100 ps, which is fast compared to the inverse tunnel cou- 
pling h/J, but still adiabatic with respect to transitions 
to higher Bloch bands. The final lattice depths were in 
the Mott-insulating regime, close to the critical point. 
After a variable evolution time, we "froze" the density 
distribution of the many-body state by rapidly raising 
the lattice depth in all directions to ~ 80 E r . Finally, the 
atoms were detected by fluorescence imaging using a mi- 
croscope objective with a resolution on the order of the 
lattice spacing and a reconstruction algorithm extracted 
the occupation number at each lattice site |19| . Because 
inelastic light-assisted collisions during the imaging lead 
to a rapid loss of atom pairs, we directly detected the 
parity of the occupation number. 

Our experimental results for the time evolution of the 
two-point parity correlations after a quench to U/J = 
9.0(3) show a clear positive signal propagating with in- 
creasing time to larger distances d (Fig. |2j. In addition, 
the propagation velocity of the correlation signal is con- 
stant over the range 2 < d < 6 (inset of Fig.|2|. We found 
similar dynamics also for quenches to U/J = 5.0(2) and 
7.0(3) (Fig. [4]). We note that the observed signal can- 
not be attributed to a simple density wave because such 
an excitation would result in (sjSj+d) — (sj)(sj+d)- We 
compared the experimental results to numerical simula- 
tions of an infinite, homogeneous system at T = using 
the adaptive time-dependent density matrix renormal- 
ization group [24, 25J (£-DMRG). In the simulation, the 
initial and final interaction strengths were fixed at the 
experimentally determined values and the quench was 
considered instantaneous, at t = 0. We found remark- 
able agreement between the experiment and theory over 
all explored distances and times, despite the finite tem- 
perature and the harmonic confinement with frequency 
v = 68(1) Hz that characterise the experimental system. 
The observed dynamics is also qualitatively reproduced 
by our analytical model for U/J = 9.0. For lower val- 
ues of U/J, however, the model breaks down due to the 
increasing number of quasiparticles. 

We extracted the propagation velocity v from the time 
of the correlation peak as a function of the distance 
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FIG. 2. Time evolution of the two-point parity cor- 
relations. After the quench, a positive correlation signal 
propagates with increasing time to larger distances. The ex- 
perimental values for a quench from U/J = 40 to U / J — 9.0 
(circles) are in good agreement with the corresponding numer- 
ical simulation for an infinite, homogeneous system at zero 
temperature (continuous line). Our analytical model (dashed 
line) also qualitatively reproduces the observed dynamics. In- 
set: Experimental data displayed as a colormap, revealing the 
propagation of the correlation signal with a well defined ve- 
locity. The experimental values result from the average over 
the central N sites of more than 1000 chains, where N equals 
80% of the length of each chain. Error bars represent the 
standard deviation. 



d (Fig. [3p). A linear fit restricted to 2 < d < 6 
yields v x h/(Ja lat ) = 5.0(2), 5.6(5) and 5.0(2) for U/J = 
5.0(2), 7.0(3) and 9.0(3), respectively. The points for 
d = 1 were excluded from the fit, as they result from 
the interference between propagating and bound quasi- 
particle pairs (see Eq. 1). A comparison of the exper- 
imental velocities with the ones obtained from numer- 
ical simulations (Fig. ^p) shows agreement within the 
error bars. The measured velocities can also be com- 
pared with two limiting cases: On the one hand, they 
are significantly larger than the spreading velocity of 
non-interacting particles, v = 4 Ja\ at /h, and twice the 
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FIG. 3. Propagation velocity, a, Determination of the 
propagation velocity for the quenches to U/J = 5.0 (trian- 
gles), 7.0 (squares) and 9.0 (circles). The time of the max- 
imum of the correlation signal is obtained from fits to the 
traces Cd (t) • Error bars represent the 68 % confidence inter- 
val of these fits. We then extract the propagation velocities 
from weigthed linear fits restricted to 2 < d < 6 (lines). The 
data for U/J — 5.0 and 7.0 have been offset horizontally for 
clarity, b, Comparison of the experimental velocities (circles) 
to the ones obtained from numerical simulations for an infi- 
nite, homogeneous system at zero temperature (shaded area). 
The shaded area and the vertical error bars denote the 68 % 
confidence interval of the fit. The horizontal error bars rep- 
resent the uncertainty due to the calibration of the lattice 
depth. The black line corresponds to the bound v ma * pre- 
dicted by our effective model (the fading indicates the break 
down of this model). The arrows mark the maximum velocity 
expected in the non-interacting case (left) and the asymptotic 
value derived from our model when U/J^oo (right). 



velocity of sound in the superfluid phase [26]; on the 
other hand, they remain below the maximum velocity 
y max « (6Jai at /fi) [l — 16J 2 /(9U 2 )] predicted by our 
analytical model, that can be interpreted as a Lieb- 
Robinson bound (Fig. [3|)). In the limit U/J — » oo, this 
bound corresponds to doublons and holons propagat- 
ing with the respective group velocities 4 Ja\ at / h and 
2 Jd\at / fi- The higher velocity of doublons simply reflects 
their Bose-enhanced tunnel coupling. 

In conclusion, we have presented the first experimen- 
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tal observation of an effective light cone for the spread- 
ing of correlations in an interacting quantum many-body 
system. Although the observed dynamics can be under- 
stood within a fermionic quasiparticle picture valid deep 
in the Mott insulating regime, we note that our exper- 
imental data cover a region close to the critical point, 
for which only ab-initio numerical simulations are avail- 
able so far |8J. Our work opens interesting perspectives, 
such as revealing the entanglement carried by the quasi- 
particle pairs or investigating the quantum dynamics in 
higher dimensions, where little is known about Lieb- 
Robinson bounds and the scaling of entanglement. For 
example, the experiment can be extended to study cor- 
relation propagation in two dimensions, where existing 
numerical and analytical approaches suffer from severe 
limitations. Furthermore, the production rate of excita- 
tions and the domain formation when tuning the effective 
interaction strength slowly across the critical point can 
be investigated, thereby exploring a quantum analog to 
the Kibble-Zurek mechanism [611271125]. 
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APPENDIX 
Quenches to U/J = 5.0 and 7.0 

We also recorded the time evolution of the two-point 
parity correlations (2) after quenches to U/J = 5.0(2) 
and 7.0(3), and compared the experimental results to 
DMRG simulations of an infinite, homogeneous system 
at zero temperature (Fig. [4|. The experimental se- 
quence was identical to the one we used for the quench 
to U/J = 9.0(3), apart from the different end point of 
the quench. The data presented here are those used in 
Fig.© 



Quasiparticle model 

In the Bose-Hubbard model, bosonic atoms in an op- 
tical lattice are confined to a single Bloch band and obey 
the Hamiltonian 



number of atoms at that site. The model is entirely 
parametrised by the effective interaction strength U/J. 

In order to analytically treat the time evolution of cor- 
relations after a sudden decrease of U/J, we developed 
a novel approach based on fermionized quasiparticles. 
The initial state being close to a Fock state with one 
atom per lattice site, an effective description of the evo- 
lution at sufficiently large final interaction strengths can 
be obtained within a local basis formed by empty states, 
|o)j? sm gly occupied states, an d doubly occupied 

states, | * )■. Using generalised Jordan-Wigner transfor- 
mations [29], we then introduced fermionic creation op- 
erators for the excess particles, d)- \ ° 9 ) ■ — ^ | \ )-, and the 

holes, hj\°)j —> | ° as well as the corresponding anni- 
hilation operators. Within the truncated Hilbert space, 
the original Hamiltonian ([3| can be exactly written in 
terms of these operators: 



H = J2 ? { ~ 2 Jd] d j+1 -Jh] +1 hj 



jV2(d] h) +1 



hj dj+i 
U 



h.c 



+ +nhj)}v, (4) 



with hd j = djdj and n^j = hjhj. The complexity of the 







— i — i — i — i — | — i — i — i — | — i — i — i — p 



where dj and Oj represent the annihilation and creation 
operator of an atom at site j and hj = d)dj counts the model is hidden in the projector V = U^I - n dJ h hJ ) 

that eliminates the unphysical situation of having an ex- 
cess particle and a hole at the same site (/ is the identity 
operator). As multiple occupancies of equal species are 
naturally avoided due to their statistics, one still obtains 
a good description of the system when setting V —> I, 
provided the density of excitations (hdj(t) + hhj(t)) re- 
mains low. This is in contrast to the usual bosonic rep- 
resentations [30, 31J. 

The Hamiltonian Q with V —± I is quadratic and can 
be diagonalised by a Bogolyubov transformation. The 
eigenmodes are doublons and holons with well defined 
momentum k: 




d = 6 



7d,fc 

7h,-fc 



u(k) 4 -| 
u(k)hl k 



v(k) h- k 
- v(k) d k 



(5) 
(6) 



with 



0.4 0.8 1.2 1.6 0.4 0.8 1.2 
time t {h/J) time t (h/J) 

FIG. 4. Time evolution of the two-point parity corre- 
lations. Left panel: quench to U/J = 5.0(2). Right panel: 
quench to U/J = 7.0(3). The circles indicate the correlations 
measured experimentally and the line is derived from the nu- 
merical simulations for an infinite, homogeneous system at 
zero temperature. The experimental and numerical values 
were obtained in the same way as described in the legend of 
Fig. [2] and in the Methods Summary section. 



u(k) = cos[0(fc)/2] , v(k) = i sin[0(fc)/2] 

and g(fc) = atan [ -^^^ 
U — 6Jcos(ka\ at ) 

Their respective eigenenergies are given by 

ed(fc) = -J cos(/cai at ) 



(7) 



1 



[U - 6Jcos(te lat )] 2 + 32J 2 sin 2 (te lat ) , (8) 



e h (k) = J cos(fcai a t) 
1 



[U - 6J cos(te lat )] 2 + 32J 2 sin^feaut) . (9) 
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Within this model, the initial state |^o) evolves in time 
according to: 

W)) = z- iiItlh IV'o) (io) 

•7l fc 7U}l^c), (11) 

with 

u(k) = u(k)u (k) - v(k)v (k) , (12) 

v(k) = v(k)u (k) - u(k)v (k) . (13) 

Here the subscript "0" denotes quantities correspond- 
ing to the initial interaction strength, whereas no la- 
bel is used for the quantities corresponding to the fi- 
nal interaction strength. Further, |vac) represents the 
quasiparticle vacuum at the final interaction strength 
(7d,fc|vac) = 7h,/c|vac) = 0). Equation (10) describes 



wave packets of entangled quasiparticle pairs that travel 
in opposite directions with different velocities. One can 
extract a maximal velocity for the spreading of the cor- 
relations from the dispersion relation of a quasiparticle 
pair (black line in Fig. ffi 



^max — 



- max 

h k 



' d_ 
dk 



ed{k)+e h {k)\ 



(14) 



Additionally, we derived from equation (10) the time 
evolution of the two-point parity correlations. In agree- 
ment with DMRG simulations, it displays a clear pos- 
itive signal, the position of which increases with time 
(Fig. [2]). We extracted the instantaneous propagation ve- 
locity v(do) from a linear fit through the signal positions 
do < d < do +4 (Fig. [5). The case do — 2 corresponds to 
the data shown in Fig. 3] At short distances, we find ve- 
locities about 10 % smaller than v max , in good agreement 
with the velocities measured experimentally. At large dis- 
tances, the velocity converges algebraically to v max . This 
latter behaviour can be understood from the expansion 
of the correlation functions to first order in J/U ', which 
can be expressed in terms of the Airy function Ai(x): 



c d (t) 



27/3^2/3^ 
3Ut 



Ai 



(2/d)^ 3 (d-6Jt/h) 



(15) 

We checked the validity of our model of freely prop- 
agating fermionic quasiparticles by comparing it with 
DMRG simulations. The propagation velocity of the two- 
point parity correlations is very accurate in the strongly 
interacting limit (e.g. U/J = 20), and remains in fairly 
good agreement down to U/J = 9, where the quasiparti- 
cle density is about 0.1 per site (Fig. |5|. At lower inter- 
action strengths, we found that the Gutzwiller approxi- 
mation breaks down. Nevertheless, the experiment and 
the simulations show that the spreading of correlations 
is still characterised by a well defined velocity, for which 
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FIG. 5. Instantaneous propagation velocity. We com- 
pare the instantaneous velocity v(do) predicted by our ana- 
lytical model (points) with the one derived from numerical 
simulations (circles). The agreement is excellent at U/J = 20 
(left panel) and qualitatively good at U/J = 9 (right panel). 
Error bars denote the 68 % confidence interval of the fit. The 
dashed line represents the asymptotic expression (15). Ar- 
rows point to the maximum velocity v ma ^ at the given inter- 
action strength. The signal position was obtained using the 
procedure described in the Methods Summary section for the 
numerical simulations. 



^max remains a relevant upper bound. We verified nu- 
merically that the truncation of the local Hilbert space 
to three states is reasonable down to U/J~6. 



Calibration of the lattice depth 

We calibrated the lattice depths by performing ampli- 
tude modulation spectroscopy of the transition between 
the zeroth and second Bloch band on a Id degenerate 
gas for the x- and 7/-axes, and on a 2d degenerate gas for 
the z-axis. We estimate the calibration uncertainty to be 
1-2%. 



Bose— Hubbard parameters 

For a given lattice depth V, we calculated the tunnel 
coupling and the on-site interaction of the Bose-Hubard 
model ([3]) in the single-particle picture using their ex- 
pressions as overlap integrals of the Wannier functions. 
In Table [T| we provide the values of V, J and U for the ef- 
fective interaction strengths mentioned in the main text. 



Ramp of the lattice depth for the quench 

The ramp of the lattice depth for the quench follows 
the functional form: 



V(t) = V + (V- V ) 



-t/r _ 1 
-T/t _ I 



(16) 



7 



U/J 


V (Er) 


J/h (s- 1 ) 


U/h (Hz) 


40(2) 


13.5(2) 


113(5) 


724(4) 


9.0(3) 


7.70(11) 


422(11) 


607(3) 


7.0(3) 


6.85(11) 


522(14) 


584(3) 


5.0(2) 


5.75(13) 


691(23) 


550(4) 



TABLE I. Lattice depth and Bose-Hubbard parameters for 
the data reported in the main text. The lattice depth in the 
perpendicular axes is 20.0(5) E v . The uncertainties on J, U 
and U / J follow from the calibration uncertainty on the lattice 
depth (see Methods Summary). 



U/J 


T (us) 


T (US) 


9.0 


100 


130 


7.0 


100 


120 


5.0 


150 


160 



TABLE II. Lattice depth ramp parameters used for the 
quenches reported in the main text, as defined in equation 
{I6f. 

Here Vb = 13.5 E T is the initial lattice depth, T is the 
duration of the ramp and r determines the rate at which 
the lattice depth is decreased. These parameters have 
to be chosen such that the ramp is fast compared to the 
time scale of the dynamics following the quench, which 
is given by h/J, but adiabatic with regard to transitions 
to higher Bloch bands. The latter condition requires the 
parameter V = (h/ A 2 )|dA/d£| to be much smaller than 1, 
where A is the energy gap between the two Bloch bands 
considered. For each quench, we have chosen the ramp 
parameters such that T is the shortest ramp duration 



compatible with T < 0.3 (Table [IT]). 

Numerical simulations showed that the origin t = 
of the time evolution can be defined as the moment 
when the effective interaction strength reaches the value 
U I J — 17. We used the same phenomenological criterion 
to locate the moment at which the dynamics stops when 
raising the lattice depth to ~ 80£? r . 



Determination of the time of the correlation peak 

We extracted the time of the maximum of the corre- 
lation signal as a function of the distance, for both the 
experiment and the theory, by fitting an offset-free gaus- 
sian profile to the traces Cd(t). For the numerical data, 
we filtered out frequency components above 3 J/h prior 
to the fit in order to isolate the envelope of the signal. For 
the experimental data, we fixed the width of the gaussian 
profile to the value obtained from the numerical data, 
keeping only the amplitude and time as free parameters. 

Numerical simulations 

The numerical simulations relied on matrix product 
states [32J (MPS) to represent infinite homogenous sys- 
tems [33| [34]. Initial states were obtained using the 
DMRG algorithm [35J. For the time evolution, we used a 
second-order Suzuki-Trotter decomposition [24[ |25| . We 
achieved quasi-exact results on the relevant time scale 
tJ/H ~ 2 by choosing a small enough Trotter time step 
(~ 0.002) and retaining a few thousand states. 



